#include "geometry.h"
class quaternion ?v ?a () =
object 

	val mutable u = new vector ()
	val mutable q = 1.0

	method getU = u
	method getQ = q

	method setU uu = u <- uu
	method setQ qq = q <- qq


	method log = let n = u#norm in
		if n < 1e-12 then u#mul (2.0 /. q) else  u#mul ((atan2 (2.0 *. n *. q) (q *. q -. n *. n )) /. n)

	method reset = u#reset ; q <- 1.0


	method mul (b:quaternion) = let v = new vector ~xx:(q *. (b#getU)#getX +. u#getX *. b#getQ +.  u#getY *. (b#getU)#getZ -. u#getZ *. (b#getU)#getY) 
		~yy:(q *. (b#getU)#getY +. u#getY *. b#getQ +.  u#getZ *. (b#getU)#getX -. u#getX *. (b#getU)#getZ)
		~zz:(q *. (b#getU)#getZ +. u#getZ *. b#getQ +.  u#getX *. (b#getU)#getY -. u#getY *. (b#getU)#getX) ()
		and qq = q *. b#getQ -. u#getX *. (b#getU)#getX -. u#getY *. (b#getU)#getY -. u#getZ *. (b#getU)#getZ 
			in let quat = new quaternion () in quat#setU v ; quat#setQ qq; quat

	method rotate (b:vector) =
		let a = u#getY *. b#getZ -. u#getZ *. b#getY +. q *. b#getX and
		b = u#getZ *. b#getX -. u#getX *. b#getZ +. q *. b#getY and
		c = u#getX *. b#getY -. u#getY *. b#getX +. q *. b#getZ and
		d = -. u#getX *. b#getX -. u#getY *. b#getY +. u#getZ *. b#getZ
		in new vector ~xx:(-. u#getZ *. b +. u#getY *. c -. u#getX *. d +. q *. a)
					~yy:(-. u#getX *. c +. u#getZ *. a -. u#getY *. d +. q *. b)
					~zz:(-. u#getY *. a +. u#getX *. b -. u#getZ *. d +. q *. c) ()

	initializer
		match (v, a) with
		(None, None) -> ()
		| (Some v, None) -> u <- v
		| (None, Some a) -> q <- a
		| (Some v, Some a) -> let halfa = 0.5 *. a in let sa = sin halfa  in u <- v#mul sa; q <- cos halfa

end;;

type triangle = {a:vector;b:vector;c:vector;n:vector};;
type mesh = {tList: triangle list; pList: vector list};;
